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Abstract 

An  efficient  full-field  method  of  computing  the  local  and  homogenized 
macroscopic  responses  of  elasto-plastic  polycrystalline  microstructnres  based 
on  the  fast  Fourier  transform  (FFT)  algorithm  is  presented.  This  approach 
takes  realistic  microstructure  images  as  the  input  and  estimates  the  mechani¬ 
cal  response/properties  of  polycrystal  microstructnres  under  periodic  bound¬ 
ary  conditions  without  requiring  complex  discretization.  Effective  stress- 
strain  response,  local  mechanical  response  fields  and  crystallographic  texture 
of  deformed  microstructnres  are  examined.  Special  interest  is  given  to  the 
fatigue  indicator  parameters  (FIPs)  of  INIOO  (nickel-based  superalloy).  This 
approach  accounts  for  both  intergranular  and  intragranular  interactions  of  an 
elasto-plastic  heterogeneous  medium  and  therefore  provides  accurate  predic¬ 
tion  of  the  mechanical  response  and  texture  evolution.  The  elastic  and  plastic 
responses  are  computed  separately  by  satisfying  two  forms  of  the  equilibrium 
equations.  A  multi-grid  strategy  is  also  adopted  to  capture  the  heterogeneous 
deformation  of  microstructnres.  The  obtained  results  are  compared  with  a 
widely  used  crystal  plasticity  finite  element  approach.  The  gained  compu¬ 
tational  efficiency  of  full-field  polycrystalline  microstructure  simulation  is 
prominent. 
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1.  Introduction 

Numerical  prediction  of  effective  and  local  mechanical  behavior  of  poly- 
crystalline  materials  has  received  great  attention  in  the  computational  ma¬ 
terials  science  community  starting  with  the  pioneering  work  on  crystal  plas¬ 
ticity  using  physically-based  models  by  Sachs  [1]  and  Taylor  [2],  The  Sachs 
model,  assumes  that  each  grain  in  the  polycrystalline  aggregate  is  subjected 
to  the  same  stress  and  provides  lower  bound  prediction  to  the  effective  stress- 
strain  response.  This  model  was  extended  to  deal  with  large  visco-plastic  and 
elasto- viscoplastic  deformations  and  is  known  as  the  lower  bound  model  [3]. 
On  the  other  hand,  the  Taylor  model  also  known  as  the  upper  bound  model 
assumes  uniform  strain  in  the  polycrystal  (equal  to  the  macroscopic  value), 
and  provides  rigid  prediction  to  the  effective  stress-strain  response  of  the 
microstructure.  Extensions  of  the  Taylor  model  to  elasto-plastic  [4],  visco¬ 
plastic  [5],  and  hnite  elasto- viscoplastic  deformations  [6]  were  proposed.  Al¬ 
though  not  accounting  for  the  interactions  inside  the  microstructure,  the 
Taylor  model  has  been  widely  used  in  computational  crystal  plasticity  com¬ 
munity  for  its  simplicity  and  high  computational  efficiency  [7,  8,  9,  10]. 

The  above  models,  however,  only  give  approximate  (or  more  precisely, 
extreme)  estimation  of  the  effective  mechanical  behavior  of  polycrystalline 
microstructures  by  assuming  homogeneous  helds  over  the  domain.  The  local 
mechanical  response  of  deformed  microstructures  cannot  be  accurately  cap¬ 
tured.  To  account  for  intergranular  interaction  during  deformation,  a  more 
sophisticated  and  popular  homogenization  strategy,  self-consistent  method, 
has  been  developed.  The  formulation  of  self-consistent  models  is  based  on 
the  solution  of  the  problem  of  an  ellipsoidal  inclusion  embedded  in  an  inhnite 
homogeneous  equivalent  medium.  The  inclusion  is  taken  to  be  an  individual 
grain  while  the  homogeneous  medium  represents  the  equivalent  polycrys- 
talline  aggregate.  Each  inclusion  (or  grain)  is  taken  as  an  averaged  medium, 
the  heterogeneity  within  which  is  not  considered.  The  hrst  attempt  to  model 
the  overall  elasto-plastic  behavior  of  polycrystals  through  the  self-consistent 
approach  was  proposed  by  Kroner  [11]  based  on  the  use  of  Eshelby’s  solu¬ 
tion  [12].  This  model  neglects  plastic  interactions  between  the  inclusion  and 
the  surrounding  matrix.  An  improvement  was  introduced  by  Hill  [13]  to  ac¬ 
count  for  the  plastic  interaction  using  an  incremental  formulation  based  on 
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the  linearization  of  the  local  constitutive  equations.  This  model  was  later  im¬ 
plemented  for  polycrystals  [14],  Hill’s  incremental  approach,  which  is  known 
as  the  secant  model,  was  extended  to  study  polycrystals  at  large  elastic- 
plastic  deformations  in  [15].  A  non- incremental  self-consistent  approach  by 
using  the  tangent  modulus-based  formulation  was  proposed  in  [16].  This  work 
has  been  used  to  simulate  large  visco-plastic  deformations  in  various  crys¬ 
talline  materials  with  full  anisotropic  overall  tangent  modulus  [17,  18].  These 
works  were  further  extended  to  a  general  elasto- viscoplastic  model  in  [19].  An 
“affine”  model,  which  yields  softer  predictions,  was  later  proposed  [20,  21] 
to  study  rate-dependent  elasto-plastic  behavior  of  polycrystals.  All  these 
developments  are  based  on  linearization  schemes  that,  at  grain  level,  only 
make  use  of  information  on  held  averages,  disregarding  higher-order  statisti¬ 
cal  information  inside  the  grain.  Improvements  taking  consideration  of  the 
second-order  moment  of  the  held  huctuation  in  grains  were  also  developed 
in  [22,  23,  24]  to  account  for  intragrain  heterogeneities. 

Self-consistent  models  are  widely  adopted  to  study  the  ehective  mechan¬ 
ical  response  and  texture  evolution  of  polycrystalline  microstructures.  Al¬ 
though  they  give  more  reliable  predictions  to  macroscopic  properties  than 
the  Sachs  and  Taylor  models,  these  “mean-held”  methods  are  still  not  capa¬ 
ble  of  estimating  local  micromechanical  helds  appropriately  since  a  homoge¬ 
neous  grain  is  assumed.  To  this  end,  full-held  simulations  that  interrogate 
poly  crystalline  microstructures  with  intracrystalline  resolution  are  developed 
with  the  evolution  of  computing  power.  The  most  popular  full-held  model 
for  polycrystals  nowadays  is  the  crystal  plasticity  hnite  element  method,  in 
which  a  variational  solution  is  achieved  for  the  force  equilibrium  and  displace¬ 
ment  compatibility  using  the  principle  of  virtual  work  for  a  domain  that  is 
discretized  into  hnite  elements.  Several  crystal  plasticity  hnite  element  mod¬ 
els  have  been  developed  incorporating  diherent  physics-based  constitutive 
laws  since  the  early  work  by  Peirce  et  al.  [25].  For  a  detailed  overview,  the 
reader  is  referred  to  [26],  where  theories  and  applications  of  crystal  plasticity 
hnite  element  models  are  reviewed.  In  our  earlier  work,  a  rate-independent 
constitutive  model  [27]  along  with  the  grain  size  ehect  described  by  lattice 
incompatibility  [28]  were  implemented  to  study  the  mechanical  response  of 
poly  crystalline  FCC  microstructures  discretized  by  conforming  grids  [29].  A 
rate-dependent  nickel-based  superalloy  constitutive  model  [30]  was  later  em¬ 
bedded  in  our  crystal  plasticity  hnite  element  framework  to  study  fatigue 
properties  and  their  variability  induced  by  microstructure  uncertainties  [31]. 

The  difficulty  in  hnite  element  meshing  coupled  with  the  large  number  of 
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degrees  of  freedom,  however,  limits  the  size  of  the  microstructures  that  can  be 
investigated  by  hnite  element  simulations.  In  addition,  the  high  computation 
cost  of  crystal  plasticity  hnite  element  analysis  also  prevents  it  from  being 
used  as  the  deterministic  solver  in  stochastic  simulations  or  as  a  point  simu¬ 
lator  in  multiscale  simulations.  Another  efficient  full-held  approach  based  on 
Green’s  functions  and  fast  Fourier  transform  (FFT)  has  been  proposed  as  an 
alternative  to  the  hnite  element  method  for  solving  the  governing  equations 
for  periodic  heterogeneous  media  [32,  33].  This  fast  Fourier  transform-based 
method  was  originally  proposed  as  a  fast  algorithm  to  compute  the  elastic 
and  elasto-plastic  ehective  and  local  response  of  composites  with  isotropic 
components  [34,  35,  36].  This  approach  was  also  adapted  to  deal  with  poly- 
crystalline  microstructures  using  a  visco-plastic  constitutive  model  [32].  Lo¬ 
cal  and  macro  mechanical  responses,  as  well  as  the  texture  evolution,  of  2D 
and  3D  realistic  polycrystals  were  stndied  in  a  series  of  works  [32,  37,  38]. 
Comparison  with  the  self-consistent  method  [39]  and  hnite  element  simula¬ 
tions  [40,  41]  were  conducted  and  showed  the  accnracy  and  efficiency  of  the 
fast  Fonrier  transform  scheme.  Recent  attempts  to  conple  the  fast  Fourier 
transform-based  model  with  hnite  elements  is  presented  in  [42].  The  major 
merit  of  the  crystal  plasticity  fast  Fourier  transform-based  method  lies  in  the 
following  aspects: 

1.  It  directly  works  on  pixelized  microstructure  images  without  requiring 
sophisticated  discretization. 

2.  It  is  a  full-held  method  that  accurately  investigates  the  global  and 
local  mechanical  behavior  of  microstructures  by  accounting  for  both 
intergrain  and  intragrain  interactions. 

3.  It  has  better  nnmerical  performance  than  the  crystal  plasticity  hnite 
element  method  for  the  same  spatial  resolntion  without  sacrihcing  ac¬ 
curacy  [40]. 

Recent  crystal  plasticity  fast  Fonrier  transform  models  neglect  elastic 
behavior  and  consider  only  visco-plastic  constitutive  behavior.  For  a  mi¬ 
crostructure  under  large  plastic  deformation,  this  approximation  is  valid  as 
elastic  strain  is  usnally  very  small.  In  certain  cases,  however,  both  elastic 
and  plastic  mechanisms  are  important,  e.g.  fatigue  analysis  of  turbine  engine 
components  nnder  cyclic  loading.  Driven  by  this  reqnirement,  we  here  intro- 
dnce  a  new  fast  Fourier  transform-based  crystal  plasticity  model  with  the  in¬ 
corporation  of  elasto-plastic  constitutive  relations  (we  call  it,  in  abbreviation, 
CEPFFT).  The  total  strain  rate  is  additively  decomposed  into  an  elastic  part 
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and  a  plastic  part.  The  elastic  and  plastic  responses  are  computed  separately 
using  the  fast  Fourier  transform-based  method  and  combined  to  update  the 
stress  held.  Through  a  series  of  benchmark  examples,  it  is  shown  that  the 
effective  and  local  mechanical  responses  predicted  by  CEPFFT  are  in  good 
agreement  with  the  crystal  plasticity  hnite  element  results.  A  constitutive 
model  of  INIOO  [30]  is  also  employed  to  study  the  microstructure-sensitive 
fatigue  indicator  parameters  of  such  Ni-based  superalloy  under  cyclic  load¬ 
ing.  Comparison  with  hnite  element  solutions  shows  great  consistence  of 
the  two  methods.  Besides,  we  analyzed  the  performance  of  the  fast  Fourier 
transform  based  simulator  with  pure  visco-plastic  model  implemented  in  two 
diherent  ways  (basic  formulation  and  augmented  Lagrangian  approach).  It  is 
observed  that  both  algorithms  work  equally  well  for  crystal  plasticity  prob¬ 
lems.  A  multi-grid  strategy  separating  the  computation  and  material  grids 
based  on  the  particle-in-cell  method  [43]  is  also  employed  and  the  obtained 
results  are  compared  with  those  using  a  single  grid  strategy  [40]. 

The  organization  of  the  paper  is  as  follows.  In  Section  2,  the  formulation 
of  CEPFFT  algorithm  is  introduced.  We  start  with  pure  elastic  and  visco¬ 
plastic  situations  and  extend  to  elasto-plastic  problem.  A  complete  boundary 
value  problem  along  with  the  solution  strategy  and  elasto-plastic  constitutive 
relation  is  described.  Specihc  constitutive  models  that  are  adopted  in  the 
paper  are  presented  in  Section  3.  The  input  model  and  microstructure  update 
strategy  are  described  in  Section  4.  Facts  about  the  crystal  plasticity  hnite 
element  method  are  briehy  reviewed  in  Section  5.  Numerical  examples  are 
demonstrated  in  Section  6,  where  comparisons  between  the  various  methods 
for  several  benchmark  problems  are  shown.  Error  analysis  and  computational 
efficiency  tests  are  conducted.  The  multi-grid  and  single-grid  results  are  also 
presented.  Finally,  a  brief  summary  is  provided  in  Section  7. 

2.  Crystal  Elasto-Plastic  Fast  Fourier  Transform  Simulator 

In  this  section,  we  address  the  solution  of  the  boundary  value  problem 
dehning  the  deformation  of  elasto-plastic  polycrystalline  microstructures  us¬ 
ing  Green’s  function  method,  in  which  any  point  in  the  domain  can  be  con¬ 
sidered  as  an  inclusion  embedded  in  a  homogeneous  reference  medium.  The 
local  mechanical  response  of  the  heterogeneous  medium  can  be  calculated  as 
a  convolution  integral  between  the  Green’s  function  associated  with  the  lin¬ 
ear  reference  homogeneous  medium  and  the  actual  heterogeneous  held  [33]. 
All  the  local  quantities  can  be  written  as  the  summation  of  a  mean  value 
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and  a  fluctuation  indicated  by  a  symbol.  Usually,  the  representative  vol¬ 
ume  element  of  bulk  microstructures  (not  on  the  surface)  is  modeled  to  be 
periodic  and  periodic  boundary  conditions  are  applied  to  control  its  defor¬ 
mation.  In  this  case,  the  Fourier  transform  can  be  employed  to  efficiently 
solve  the  problem  in  Fourier  space,  where  the  convolution  in  the  real  space  is 
reduced  to  simple  product.  An  iterative  scheme  is  needed  to  ensure  that  the 
solution  converges  to  the  micromechanical  responses  satisfying  equilibrium 
and  compatibility  conditions. 

For  clarity,  we  will  start  with  brief  presentations  of  the  solutions  of  pure 
elastic  and  pure  visco-plastic  problems.  Then,  the  solution  strategy  of  elasto- 
plastic  problems  will  be  introduced.  For  elasto-plastic  problems,  the  strain 
rate  is  coupled  with  stress  and  stress  rate  leading  to  challenges  in  using  the 
Green’s  function  method.  To  address  this,  we  will  introduce  a  solution  strat¬ 
egy  in  which  the  elastic  and  plastic  responses  are  computed  simultaneously 
but  separately. 

2.1.  Solution  of  Crystal  Elastic  Boundary  Value  Problems 

In  pure  crystal  elastic  problems,  the  total  strain  rate  s  is  equal  to  the 
elastic  strain  rate  and  the  stress  rate  is  linearly  related  to  the  strain  rate 
through  the  generalized  Hooke’s  lawd 

d-(x)  =  C"(x)  ;  £"(x),  (1) 

where  Rjn{yi) Rko{'^) Rip(2^)Cmnop  is  the  local  elastic  stiffness 

tensor  represented  in  the  sample  coordinate  system  that  relates  the  crystal 
lattice  frame  via  a  rotation  matrix  R(x)  determined  by  the  orientation  r(x) 
of  the  crystal  at  position  x.  C®  is  the  elastic  stiffness  modulus  in  the  lattice 
coordinate  system. 

For  a  microstructure  subjected  to  periodic  boundary  conditions  with  an 
imposed  average  velocity  gradient  L  =  VV,  the  local  equilibrium  equation, 
represented  in  terms  of  stress  rate,  needs  to  be  satisfied  at  any  point  x  within 
the  microstructure  domain  B.  The  complete  boundary  value  problem  is  de¬ 
fined  as: 

V  •  <t(x)  =  V  ■  d-(x)  =  0  Vx  G  H, 
is  periodic,  &  ■  n  is  antiperiodic  on  dB,  (2) 


^We  state  the  elastic  problem  in  terms  of  stress  and  strain  rates  for  exploring  similarities 
with  the  formulation  of  the  elasto-plastic  problem  in  Section  2.2. 
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where  v®(x)  =  v®(x)  —  L  ■  x  is  the  velocity  fluctuation  (deviation  of  the  local 
velocity  v®(x)  from  the  mean  velocity  V)  at  x  induced  by  the  microstructure 
heterogeneity.  The  superscript  e  indicates  that  the  response  stems  from 
elastic  deformation.  Our  goal  is  to  compute  the  velocities  and  their  gradients 
of  all  material  points  that  satisfy  the  above  governing  equations  and  use  them 
to  evaluate  the  strain  and  stress  responses  over  the  microstructure  domain. 
To  this  end,  we  can  write  the  local  stress  cr(x)  as  the  sum  of  two  terms: 

cr(x)  =C°:s'=(v'=(x))  +  0'=(x).  (3) 

In  Eq.  (3),  C°  is  the  stiffness  modulus  of  a  linear  homogeneous  reference 
medium,  in  which  point  x  is  imbedded.  In  this  elastic  problem,  C'^  is  selected 
as  the  averaged  elastic  modulus  C®  over  the  microstructure  domain: 

C»  =  c«  =  (C')^  =  i^C'(x)dx.  (4) 

The  second  term,  </>'^(x),  in  Eq.  (3)  is  the  periodic  polarization  field  de¬ 
fined  as 

0"(x)  =  &(x)  -  :  £"(v^(x))  =  ^(x)  -  (5) 

Substituting  Eq.  (3)  into  the  equilibrium  equation  (Eq.  (2)),  we  obtain 

^ijki^kij  +  4>ij,j  =  0  or  +  4>ij,j  =  0-  (6) 

Representing  the  strain  rate  in  terms  of  the  velocity  gradient  |  j) , 

we  obtain: 

Qjkivhj  +  =  0  or  =  0.  (7) 

The  differential  Eq.  (7)  can  be  solved  by  means  of  the  Green’s  function 

method.  Introducing  the  Green’s  function  G®(x,  x'),  the  solution  vl{x)  takes 
the  form 

=  -  [  G^^(X  -  x')Cn,n(x')dx'.  (8) 

JB 

Substituting  Eq.  (8)  into  Eq.  (7),  leads  to: 


^eO 

'^ijkl 


x)Cn,n(x')dx' 


x')dx' 


0,  (9) 
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which  can  be  rearranged  as 


[  [-C!jklGlm,lj(^  -  x')  +  SimS{x  -  x')]  0^„^„(x')dx'  =  0.  (10) 

Jb 

Taking  (p^nn  to  be  arbitrary,  we  arrive  at  the  local  equilibrium  equation  in 
the  Green’s  function  form: 

-  ^')  +  -  ^')  =  0-  (11) 

This  equation  can  be  transformed  to  Fourier  space  where  the  convolutional 
solution  in  the  real  space  (Eq.  (8))  is  represented  by  a  simple  product.  The 
equilibrium  equation  in  Fourier  space  is  then: 

=  -S.m,  (12) 

where  ^  is  a  point  (frequency)  in  the  Fourier  space.  Solving  the  linear  system 
Eq.  (12),  we  obtain  the  Green’s  function  in  Fourier  space  Gl^  to  be 

Atk  =  (13) 

Dehning 

(14) 

and  integrating  Eq.  (8)  by  parts  while  assuming  that  the  boundary  terms 
vanish  [32],  we  can  compute  the  velocity  fluctuation  as: 

l'fc(x)  =  [  -  x')Cn(x')dx'.  (15) 

Jb 

The  velocity  and  its  gradient  fluctuations^  in  Fourier  space  are 

g(^)  =  i0G'L(OC,(O,  (16) 

Cfc(^)  =  (17) 

^Note  that  the  fluctuations  can  also  be  total  velocity  and  its  gradient,  depending  on 

the  assignment  of  the  0  frequency  term. 
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After  transforming  them  back  to  the  real  space  (e.g.  v®(x)  =  FFT  ^(v  (^))), 
the  strain  rate  and  spin  rate  flnctnations  can  be  evalnated,  respectively,  by 


chLlX) 


1 


Mj- 


(18) 


The  stress  rate  can  be  npdated  according  to  the  Hooke’s  Law  (Eq.(l)).  With 
the  npdated  stress  rate  and  strain  rate,  we  can  perform  the  next  iteration 
nntil  converged  results  are  reached. 

2.2.  Solution  of  Crystal  Visco-Plastic  Boundary  Value  Problems 

If  the  deformation  is  assumed  to  be  pure  visco-plastic  (i.e.  the  elastic 
response  is  completely  neglected,  e  =  e^).  A  nonlinear  constitutive  model  is 
employed  to  link  stress  to  strain  rate  in  the  form: 


£:^(x)  =  M^(cr(x))  ;  cr(x). 


(19) 


where  M^(<t(x))  is  a  plastic  compliance  tensor  that  is  nonlinearly  dependent 
on  stress  cr(x). 

The  solution  procedure  of  a  visco-plastic  problem  is  similar  to  that  in 
Section  2.1  except  that  this  time  the  equilibrium  equation  is  written  in  terms 
of  stress  rather  than  stress  rate: 


V  ■  cr(x) 

is  periodic. 


=  V  ■  cr(x)  =  0  Vx  G  H, 

<T  ■  n  is  antiperiodic  on  dB, 


(20) 


where  the  superscript  p  indicates  plastic  deformation  induced  response.  In¬ 
compressibility  of  plastic  deformation  also  needs  to  be  satished.  In  [32],  a 
Fourier  transform-based  algorithm  was  proposed  where  the  incompressibil¬ 
ity  condition  was  satished  by  introducing  explicitly  the  constraint  =  0. 
In  the  current  work,  the  incompressibility  condition  is  accounted  for  by 
computing  the  polarization  0^(x)  with  a  strain  rate  updated  iteratively  as 


—  \tr{£^)  as  follows: 

<#.''(x)  =  <t(x)  -  O'*  :  e’’(v''(x))  =  ct(x)  -  C'”  :  F(v»(x)), 


(21) 


where  the  stiffness  modulus  of  the  linear  homogeneous  reference  medium  is 
taken  to  be  the  averaged  plastic  modulus  over  the  microstructure 
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domain: 


C-»  =  {C-),  =  ^Jc-{a(x))dx.  (22) 

The  form  of  the  plastic  modulus  or  equivalently  the  plastic  compliance 
is  determined  by  the  specihc  plastic  constitutive  model  that  is  adopted.  In 
Section  3,  we  will  present  two  constitutive  models  to  be  used  in  the  examples. 
The  plastic  problem  can  be  solved  following  the  steps  from  Eq.  (6)  to  Eq.  (18) 
simply  by  replacing  the  superscript  e  to  p  and  using  stress  instead  of  stress 
rate. 

2.3.  Solution  of  Crystal  Elasto- Plastic  Boundary  Value  Problems 

For  elasto-plastic  problems,  the  total  strain  rate  e  is  additively  decom¬ 
posed  into  an  elastic  term  and  a  plastic  part  sP  with  trfe^)  =  0: 

£(x)  =  £®(x)  £^(x).  (23) 

Following  the  constitutive  relations  given  above,  both  stress  <t  and  stress  rate 
&  are  entangled  with  strain  rate: 

£(x)  =  M®(x)  :  <t(x)  -|-  M^(cr(x))  :  cr(x),  (24) 

where  M®  =  is  the  elastic  compliance  tensor.  The  local  mechanical 

response  at  time  t  depends  on  the  entire  loading  history  of  the  specimen. 
The  current  stress  cr  is  computed  by  cr  =  -|-  d-dt,  where  cr„  is  the  stress 

at  the  previous  time  step. 

To  follow  a  Green’s  function  approach  to  the  elasto-plastic  boundary  value 
problem,  a  proper  modulus  C°,  taking  both  elasticity  and  plasticity  into  ac¬ 
count,  needs  to  be  designed  for  the  homogeneous  reference  medium  that  can 
directly  link  either  stress  or  stress  rate  to  strain  rate.  However,  to  design  such 
an  effective  modulus  is  not  trivial.  Therefore,  we  here  propose  a  formulation 
that  solves  for  the  elastic  and  plastic  velocity  fluctuations  separately  using 
the  fast  Fourier  transform-based  algorithm  thus  avoiding  the  construction  of 
the  elasto-plastic  modulus.  The  total  velocity  gradient  at  a  single  point  is 
then  computed  by  adding  the  two  fluctuations  to  the  mean  value  VV.  After 
that,  a  nonlinear  constitutive  model  is  designed  to  update  the  stress  and 
stress  rate  given  the  total  strain  rate.  We  will  this  approach  as  CEPFFT. 

The  key  of  the  CEPFFT  approach  is  that  we  solve  simultaneously  the  two 
forms  of  the  equilibrium  equations  dehned  in  Eq.  (2)  and  (20)  for  elastic  and 
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plastic  velocity  (and  their  gradient)  fluctuations.  The  total  velocity  gradient 
can  then  be  obtained  as 


Vv(x)  =  VV  +  Vv"(x)  +  Vv^’(x),  (25) 

from  which,  the  total  strain  rate  e  can  be  calculated  as  the  symmetric  part 
of  Vv.  However,  the  portion  of  elastic  strain  rate  e®  and  plastic  strain  rate 
in  e  is  not  known.  As  a  result,  the  stress  and  stress  rate  corresponding  to 
a  given  total  strain  rate  cannot  be  directly  computed. 

An  iterative  scheme  is  designed  to  linearize  the  nonlinear  relations  among 
stress,  stress  rate,  and  strain  rates,  in  order  to  provide  a  way  of  finding  the 
elastic  and  plastic  strain  rates  given  the  total  strain  rate.  We  first  rewrite 
Eq.  (23)  as 

F  =  £'^(x)  +  £^(x)  —  £(x)  =  0.  (26) 


We  aim  at  solving  this  equation  for  the  elastic  strain  rate  with  known  e 
using  the  Newton-Raphson  scheme.  To  this  end,  Eq.  (26)  can  be  linearized 
as  follows: 


dF 


,e(i+l) 


(27) 


According  to  the  elastic  and  plastic  constitutive  models,  as  well  as  the  stress 
incremental  equation  <t  =  cr„  +  At&,  we  can  write  the  following  relations: 


deP 

dcr 

dcr 

d& 

d& 

dP 


Mr, 

=  Afll 

d& 

C", 


(28) 


where  II  is  the  fourth-order  identity  tensor  and  the  tangent  plastic  compli¬ 
ance  Mf  is  specified  by  the  particular  plastic  constitutive  model  used  (see 
Section  3.1).  Using  these  relations,  we  can  simplify  Eq.  (27)  as: 


p(i+l)('^e(i+l)^ 


Fh)(£db)  + 


Fh)(sdb)  + 


/ds" 
(dF  + 

dir" 

-TTV  + 


do-  dd-dsV  ■  ^  > 


Fh)(£db)  +  (II  +  AtMf  :  C")  :  .  (29) 
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Setting  =  0,  the  elastic  strain  rate  at  iteration  i  +  1  is  computed 


by 


^e(i+l)  ^  ^e{i)  _  ^  .  Qeyl  .  p(*)^ 

With  the  elastic  strain  rate  £®h+i)  (i+l)th  iteration  computed  from 

the  equation  above,  the  stress  rate  can  be  obtained  using  Eq.  (1). 

Therefore,  the  stress  is  updated  as  crh+i)  _  _l_  ^{i+i)  ^  with  which  the 

plastic  strain  rate  e^’h+i)  gg^j^  ^g  computed  using  the  plastic  constitutive  re¬ 
lation  Eq.  (19).  The  impressibility  is  enforced  by  setting  ^  _ 

Iteratively  updating  the  above  equations,  the  hnal  elastic  as 
well  as  the  plastic  strain  rates  can  be  computed  until  convergence  is  achieved 
(i.e.  when  F  in  Eq.  (26)  approaches  0).  The  converged  stress  rate  and  stress 
can  be  subsequently  computed  according  to  the  constitutive  model.  After 
that,  we  can  construct  the  elastic  and  plastic  polarization  helds,  respectively, 
following  Eqs.  (5)  and  (21).  By  transforming  them  to  Fourier  space,  the 
fluctuations  of  velocity  gradients  induced  by  elastic  and  plastic  deformation 
can  be  updated  using  Green’s  functions.  Inversely  transforming  these  fluc¬ 
tuations  to  real  space,  a  new  strain  rate  held  as  well  as  stress  and  stress  rate 
can  be  obtained.  The  algorithm  can  then  proceed  to  the  next  iteration.  The 
overall  CEPFFT  algorithm  is  summarized  next. 

2.4-  CEPFFT  Algorithm 

We  adopt  a  basic  fast  Fourier  transform-based  algorithm  for  implement¬ 
ing  the  CEPFFT  simulator.  The  algorithm  was  hrstly  proposed  in  [34,  35] 
to  deal  with  elastic  and  elasto-plastic  heterogeneous  composites  with  mod¬ 
erate  contrast.  It  was  also  adapted  by  Lebensohn  [32]  to  solve  elastic  and 
visco-plastic  polycrystalline  problems.  This  method  is  based  on  the  exact 
expression  of  Green’s  function  for  linear  elastic,  homogeneous  reference  ma¬ 
terial. 

Algorithm: 

1.  At  the  hrst  iteration,  start  with  an  initial  guess  of  the  total  velocity 
gradient  held  at  time  step  n  +  1-.  °Vv„+i(x)  =  Vv„(x),  Vx  G  B,  from 
which  the  local  strain  rate  can  be  computed  £(x)  =  s|/m(Vv(x)).  Then 
evaluate  the  elastic  portion  £'^(x)  and  plastic  portion  £^(x)  of  the  strain 
rate  using  the  local  constitutive  relations.  At  the  same  time,  compute 
the  initial  stress  °cr(x)  and  stress  rate  °d’(x). 
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2.  Compute  the  elastic  and  plastic  polarization  fields,  and 

respectively,  for  the  i  iteration  given  the  stress,  stress  rate  and  strain 
rate  fields: 


>"(x)  =*  d-(x)  -  £"(x), 

*0P(x)  ='  cr(x)  -  £P(x).  (31) 

3.  Transform  the  polarizations  to  Fourier  space  via  fast  Fourier  transform: 

=  FFT{^cj)%^))  and  =  F FT  cff  {-k)) . 

4.  Compute  the  velocity  gradient  fluctuation  fields  induced  by  elastic  and 
plastic  deformation,  respectively,  in  the  Fourier  space  for  the  {i  +  l)-th 
iteration. 

=  T{^)  ^  0;  and  =  0, 

=  and  '+Vv'’(0)  =  0.  (32) 

5.  Transform  the  current  velocity  gradient  fluctuation  fields  back  to  the 

real  space  through  inverse  fast  Fourier  transform,  i.e.  *’''^Vv®(x)  = 
FFT-^  and  *+iVvP(x)  =  (*+^Vv^(|)). 

6.  Compute  the  total  strain  rate  =  E  k  and  then  the 

stress  *’''^<t(x)  and  stress  rate  *+^d’(x)  fields  according  to  the  constitu¬ 
tive  model. 

7.  Check  the  error  (equilibrium  condition): 

_  (iiv  ■  ■+vir^)’-''^  _ 

||*+icr||  ||*+^d-(0)|| 

If  e  is  smaller  than  a  prescribed  error  tolerance,  the  iteration  process 
stops.  Otherwise,  return  to  step  (2)  and  proceed  to  the  next  iteration. 


Upon  convergence,  the  grain  orientations  and  the  positions  of  the  pixel 
points  are  updated  according  to  the  spin  and  velocity  gradient  fields,  respec¬ 
tively  (see  Section  4.3  for  details). 

Remark  1:  The  error  of  the  equilibrium  condition  (Step  (7))  is  mostly 
determined  by  the  resolution  of  the  image  (as  will  be  shown  later  in  the 
examples).  For  images  with  coarse  resolution,  the  error  may  converge,  with 
fluctuations,  to  a  value  larger  than  0.  This  error  is  inherently  associated 
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with  the  FFT-based  methodology.  Therefore,  a  practical  way  to  check  the 
convergence  is  to  examine  the  relevant  difference  between  the  error  at  the 
current  and  last  iterations: 


i+ig  _i  g 


Remark  2:  It  is  clear  that  the  computational  cost  of  one  iteration  step  of  the 
algorithm  discussed  above  is  doubled  that  corresponding  to  either  the  pure 
elastic  or  visco-plastic  algorithms.  If  a  proper  elasto-plastic  modulus  can  be 
designed,  the  problem  can  be  solved  using  only  one  set  of  equations,  which 
may  lead  to  higher  computational  efficiency.  A  feasible  approach  arises  by 
representing  the  local  stress  cr(x)  in  terms  of  the  total  strain  rate  £(v(x)): 

cr(x)  =  C°  ;  £(v(x))  +  0(x),  (35) 


where  C°  =  is  computed  as  the  volume  average  of  an  elasto-plastic 

modulus  derived  as  follows: 


cep  ^ 


dcr 

de 


dcr  d£®  dcr  de^ 
de®  de  de^  de 


AtC"  +  QP. 


(36) 


The  plastic  modulus  is  chosen  to  be  the  secant  modulus  that  depends 
on  the  specihcs  of  the  constitutive  model  adopted.  With  the  construction 
of  the  elasto-plastic  modulus,  the  problem  can  be  solved  following  the  same 
procedure  as  for  the  visco-plastic  problem.  The  convergence  rate  of  this 
particular  algorithm  will  be  shown  in  Section  6  that  is  slower  than  the  main 
algorithm  presented  earlier.  A  more  sophisticated  design  of  the  elasto-plastic 
modulus  is  of  great  interest. 


3.  Constitutive  Model 

The  fast  Fourier  transform  based  full  field  algorithm  can  be  adapted  to 
various  materials  assuming  appropriate  constitutive  models  are  provided.  In 
this  section,  we  will  describe  two  plastic  constitutive  models  that  are  adopted 
in  the  numerical  examples  for  single  phase  FCC  crystal  and  INIOO  superalloy. 
The  plastic  modulus  that  is  needed  in  the  CEPFFT  method  is  derived.  The 
fatigue  indicator  parameters  that  will  be  employed  to  measure  the  fatigue 
properties  of  Ni-based  superalloys  are  dehned  in  [44]. 
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3.1.  Plastic  Constitutive  Relations 

Various  plastic  constitutive  relations  have  been  developed  based  on  differ¬ 
ent  flow  rules  and  hardening  laws.  Usually,  the  nonlinear  connection  between 
plastic  strain  rate  and  Cauchy  stress,  £^(<t(x)),  is  assumed,  where  the  plas¬ 
tic  strain  rate  is  defined  as  the  sum  of  the  contributions  of  shear  rates, 
7^"^(x),  on  all  active  slip  systems: 


Ns 

sP(x)  =  ^m(“)(x)7(“)(x).  (37) 

a 

In  the  above  equation,  Ng  is  the  number  of  active  slip  systems,  and 
denotes  the  symmetric  Schmid  tensor  of  slip  system  a: 

^  (s(“)  (g)  0  ,  (38) 

where  and  are  the  slip  direction  and  slip  plane  normal  of  the  system 
a,  respectively. 

The  shear  rate  7*^“)(x)  is  determined  by  the  resolved  shear  stress  = 
;  cr  =  ;  cr  and  slip  resistance  on  slip  system  a  through  a  specihc 

flow  rule.  The  difference  between  the  various  plastic  constitutive  models  lies 
in  the  choice  of  the  flow  rule  and  hardening  law  that  dehnes  the  evolution  of 
the  slip  resistances. 

3.1.1.  A  typical  flow  rule  and  hardening  law  for  single  phase  crystalline  ma¬ 
terials 

The  hrst  model  we  introduce  here  is  a  typical  one  for  single  phase  poly¬ 
crystals.  Asaro  and  Needleman  [6]  proposed  a  widely  used  rate-dependent 
crystal  plasticity  flow  rule  for  evaluating  shearing  rate,  7^“\  on  slip  system 
a: 


(l/m) 

sign(r(“)),  (39) 

where  70  is  a  reference  rate  of  shearing  and  m  characterizes  the  material  rate 
sensitivity.  is  the  slip  resistance  of  system  a.  The  slip  resistance  evolves 
following  certain  hardening  law  that  is  calibrated  from  experimental  data  to 


r 


(a) 
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describe  the  strain  hardening  effect.  Currently,  we  adopt  a  Voce  type  model 
given  in  [40],  where  the  hardening  rate  is  computed  by 

dr  ^ 


with  the  hardening  function 


—  ^0  ff  (^1  ff  ^ih) 


0or\l 

1  —  exp 

L  V 

Ki  ) 

(41) 


where  is  a  hardening  matrix  whose  diagonal  elements  denote  self-hardening 
and  off-diagonal  elements  denote  latent  hardening.  In  the  current  work,  we 
assume  both  of  the  two  are  identical  and  equal  to  1.  kq,  ki,  6*o,  6i  are  material 
dependent  parameters.  The  cumulative  shear  T  in  a  grain  is  defined  as 


r  = 


(42) 


Following  the  above  equations,  we  can  derive  an  explicit  form  of  the 
plastic  constitutive  law  as 


e''(x)  =  M;((t(x))  :  (t(x), 


(43) 


where  is  called  the  secant  plastic  compliance  having  the  form  given  below: 


Ns 


(<^(x))  =70$^ 


m 


(«)i 


X 


m 


(“)( 


K 


(“)( 


m 


(“) 


X  ;  cr  X 


K 


(l/m-l) 


(44) 


The  corresponding  tangent  plastic  compliance  in  Eq.  (30)  at  this  spe¬ 
cific  case  is  Mf  =  The  plastic  modulus  in  Eq.  (36)  for  the  method 

highlighted  in  Remark  2  is  Given  a  plastic  strain  rate,  the 

corresponding  stress  needs  to  be  computed  following  a  iterative  scheme  (e.g. 
Newton-Raphson  method)  because  of  the  nonlinear  nature  of  the  constitutive 
model. 


3.1.2.  A  specific  flow  rule  and  hardening  law  for  two-phase  nickel-based  su¬ 
peralloy  IN  100 

In  this  work,  we  are  also  interested  in  the  fatigue  properties  of  INIOO,  a 
nickel-based  superalloy,  at  high  temperature  (650°C).  Since  the  microstruc¬ 
ture  of  superalloys  is  complex,  more  sophisticated  flow  rule  and  hardening 
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laws  are  developed  to  describe  the  micromechanical  behavior  of  two-phase 
microstructure.  In  this  work,  we  employ  the  homogeneous  constitutive  model 
developed  in  [30].  The  second  phase  (7'  precipitates)  conhguration  is  not  ex¬ 
plicitly  modeled.  Instead,  the  effects  of  the  second  phase  are  taken  into 
account  through  particular  parameters.  Cube  slip  (110){100}  systems  are 
activated  in  addition  to  octahedral  slip  (110) {111}  systems  to  take  cross  slip 
mechanism  at  high  temperatures  into  consideration.  The  rate  dependent  flow 
rule  which  estimates  the  shearing  rate  on  each  slip  system  includes  a  back 
force  term  y  for  the  modeling  of  the  Baushinger  effect  arising  principally  from 
matrix  dislocation  interaction  with  7'  phase.  The  effect  of  volume  and  size 
of  7'  precipitates  on  material  strength  is  taken  into  account  by  constitutive 
parameters.  The  constitutive  equations  are  summarized  below  and  detailed 
in  [30,  44,  45]. 

The  flow  rule  of  slip  system  a  is 
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(a) 


+ 


•  (“) 
72 


(“)  I 

I 


K 


(a) 


ni 


(a) 


D 


(a) 


n2- 


sgn(r 


(«) 


(45) 


where  7^“^  and  72°"^  are  constants  related  to  the  initial  shearing  rate  and 
is  the  drag  stress  assumed  to  be  constant.  A  ={oct,  cub}  refers  to  the 
octahedral  and  cube  slip  systems,  respectively.  The  function  {x)  returns  x  if 
a;  >  0  and  returns  0,  otherwise. 

The  evolution  of  the  slip  resistance  follows  the  Taylor  strain  harden¬ 
ing  law  determined  by  dislocation  density 


4“^  =  4“a  +  (46) 

where  at  =  0.0385  and  firnix  =  (/pi  +  fp2  +  fp3)f^Y  +  /Uy  and  are 

shear  moduli  for  7'  precipitates  and  7  matrix,  respectively.  The  magnitude 
of  Burgers  vector  is  6  =  (/pi  -h  fp2  +  fp3)by  +  fj)^.  /pi,  /p2,  fp3  are  volume 
fractions  of  primary,  secondary,  and  tertiary  7'  precipitates,  respectively,  and 
fm  =  1  — /pi  — /p2  — /p3  is  the  volume  fraction  of  7  matrix  phase, 

fv2  =  t  and  /„3  =  f  .  For  different  slip  systems,  the  initial  slip 
resistance  can  be  evaluated  by 

+  if  pi  +  /p2)t1s  ) 


(a) 

0,oct 


(h 


(a) 


0,oct/ 
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(47) 


where 


K, 


(a) 

O^cub 


+  i> 


rik 

cub 


W 


and 


+  Cps^wf^^ds  + 


APB 


^APB-ref' 


(48) 


=  hpe^ie^  +  hcb\r^h^  \  +  kseT, 


-(“) 


(49) 


in  which  T apb  is  the  anti-phase  boundary  energy  density  here  taken  be  equal 
to  TAPB-ref,  di,  i  =  1,2,3  are  the  sizes  of  precipitates,  and  dgr  is  the  grain 
size.  The  non-Schmid  stress  is  related  to  the  resolved  shear  stresses  on 
the  primary  Tpe\  cube  and  secondary  slip  systems,  respectively,  by 
parameter  hpe,  hcb,  and  hse-  Currently,  we  assume  this  term  to  be  0. 

The  dislocation  density  evolution  has  the  following  form: 


ho 


i^Zo  +  ki^x 


Zo  = 


h 


(50) 


The  evolution  of  the  back  stress 
density  and  shear  rate: 


is  also  dependent  on  dislocation 


Xa“^  =  C'x  ^V\f^mixb\f^sgn{T^°^^  - 

Vo,xZo 


(51) 


Vx  = 


Zq  +  ki^x\l Px  ^ 


where  =  123.93  -  433.98f^2  +  384.06/;|. 

Parameters  for  superalloys  at  650°C  listed  in  [30]  are  adopted.  For  addi¬ 
tional  information  about  the  constitutive  model,  the  reader  is  referred  to  [44]. 
Given  the  stress,  the  plastic  strain  rate  needs  to  be  computed  following  the 
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nonlinear  relationship  and  vice  versa.  Eq.  (30)  requires  the  evaluation  of 
tangent  plastic  compliance.  In  the  case  of  INIOO  superalloy,  the  compliance 
can  be  derived  as: 
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The  secant  plastic  modulus  required  by  Eqs.  (21)  and  (36)  is  taken  to  be 
CP  =  Mf-i  with 
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With  these  relations,  the  algorithm  proposed  in  Section  2  can  be  applied  to 
INIOO  superalloy. 


4.  Microstructure  Model 

The  solution  strategy  of  computing  the  deformation  of  polycrystalline 
microstructures  under  periodic  boundary  conditions  was  discussed  in  Sec¬ 
tion  2.  The  main  procedure  is:  (1)  compute  polarization  held;  (2)  transform 
the  polarization  held  to  Fourier  space  using  fast  Fourier  transform;  (3)  up¬ 
date  velocity  gradient  in  the  Fourier  space  and  transform  it  back  to  the  real 
space;  (4)  update  real-space  helds  (e.g.  strain  rate,  spin  rate,  stress  rate, 
etc.)  accordingly.  An  iterative  scheme  is  adopted  to  obtain  convergence. 
In  this  section,  we  would  introduce  the  digital  microstructure  model  that  is 
used  as  the  input  to  FFT-based  simulations.  The  update  strategy  of  the 
microstructure  and  crystal  orientation  during  deformation  is  also  described. 
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4.1.  Discretization 

The  input  to  the  FFT-based  (including  pure  elastic,  pure  visco-plastic, 
and  elasto-plastic)  simulators  is  pixelized  images  with  orientation  parameters 
associated  with  each  pixel  (or  voxel  for  3D).  The  pixels  or  voxels  are  the 
discretization  of  the  input  image,  which  essentially  requires  no  effort. 

To  implement  the  FFT-based  algorithm  and  solve  the  underlying  bound¬ 
ary  value  problem,  the  microstructure  is  discretized  by  a  regular  grid  con¬ 
sisting  of  Ni  X  N2  pixels  (2D  problem)  or  iVi  x  iV2  x  voxels  (3D  problem). 
Denote  Li  to  be  the  period  (edge  length)  of  the  microstructure  in  the  ith 
direction  {i  =  1,2  for  2D  and  i  =  1,2,3  for  3D).  The  coordinates  of  the 
point,  which  is  pixel  in  2D  and  voxel  in  3D,  in  the  ith  direction  is  therefore: 


(54) 


In  Fig.  1,  examples  of  2D  and  3D  grids  of  polycrystalline  microstructures 
are  shown  in  comparison  with  their  image  views.  Different  colors  in  the  mi¬ 
crostructure  represent  grains  with  distinct  orientations.  The  grain  structures 
are  generated  following  a  Voronoi  tessellation  scheme,  where  the  positions  of 
centroids  are  adjusted  to  minimize  the  interaction  forces  [46]  (see  Section  4.2). 

The  regular  discretization  grid  of  the  microstructure  determines  a  regular 
reciprocal  grid  in  the  Fourier  space,  which  makes  the  fast  Fourier  transform 
convenient.  The  i-th  direction  coordinates  of  the  points  in  the  reciprocal 
grid,  namely  frequencies,  are 


(55) 


where  i  =  1,  2  for  2D  and  i  =  1,  2,  3  for  3D.  In  the  current  work,  the  number 
of  points  in  each  dimension  is  selected  to  be  a  power  of  2  in  order  to  facilitate 
the  fast  Fourier  transform  that  is  conducted  using  the  FFTW  libraries  [47]. 


4-2.  Micro  structure  Generation 

The  microstructure  image  can  be  generated  by  several  experimental  or 
numerical  ways,  e.g.  electron  backscattered  diffraction  (EBSD),  phase  held 
method,  Monte  Carlo  grain  growth,  or  Voronoi  tessellation.  In  the  current 
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Figure  1:  (a)  The  image  representation  of  a  2D  polycrystalline  microstructure  containing 
10  grains,  (b)  The  pixel  grid  of  the  10-grain  2D  microstructure.  The  microstructure 
is  discretized  by  16  x  16  pixels,  (c)  The  image  representation  of  a  3D  polycrystalline 
microstructure  containing  64  grains,  (b)  The  voxel  grid  of  the  64-grain  3D  microstructure. 
The  microstructure  is  discretized  by  16  x  16  x  16  voxels. 


work,  we  adopt  the  Voronoi  tessellation  scheme  in  combination  with  a  “re¬ 
laxation”  strategy  to  generate  polycrystalline  microstructures  with  approxi¬ 
mately  controlled  grain  sizes.  The  algorithm  is  summarized  as  follows: 

•  Given  a  set  of  target  sizes  of  the  microstructure  constituent  grains,  we 
hrst  randomly  place  Voronoi  cell  centroids  within  the  microstructure 
domain.  The  number  of  centroid  points,  Np,  is  decided  by  Np  =  V/ fid, 
where  V  is  the  microstructure  volume  (for  3D)  or  area  (for  2D),  and 
fid  is  the  mean  volume  or  area  of  grains. 

•  A  force  function  F{i,j)  is  dehned  for  the  adjustment  of  grain  centroid 
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locations: 


F{i,j)  =  n{j,i)  *  {r{i)+r{j)  -dist{i,j)),  (56) 

where  the  force  exerted  by  point  j  on  point  i  in  the  direction 

which  is  along  the  line  starting  from  j  and  pointing  to  i.  r{i) 
is  the  radius  (assuming  spherical  (for  3D)  or  circular  (for  2D)  grains) 
of  grain  i.  dist{i,j)  is  the  distance  between  point  i  and  point  j,  which 
can  be  periodic  if  the  microstructure  is  periodic.  Function  (x)  =  x,  if 
a;  >  0,  and  (x)  =  0,  otherwise. 

•  Relaxation  step:  move  points  according  to  the  velocity  function  dehned 
as 


v(i)  =  F(*)/l/(*),  (57) 

where  F(i)  is  the  net  force  of  point  i:  F(i)  =  and  V {i)  is  the 

volume  or  area  of  the  corresponding  grain.  The  location  of  centroid  i 
is  updated  to  Pn+i(i)  =  Pn(*)  +  v(i)At,  where  At  adaptively  changes 
for  the  convenience  of  convergence. 

•  After  all  centroid  locations  are  stable  (force  balance  is  reached),  a 
Voronoi  tessellation  is  computed  using  them.  The  polycrystalline  mi¬ 
crostructure  is  therefore  constructed. 

A  schematic  description  of  relaxation  of  the  Voronoi  cell  centroids  is  de¬ 
picted  in  Fig.  2. 

4-3.  Grid  and  Texture  Update 

The  grid  after  deformation  may  become  irregular  as  material  points  move 
according  to  local  velocities.  The  exact  new  position  of  material  point  X  is 

x(X)  =  X+ (L  + Vv(X))XAf,  (58) 

where  L  =  VV  is  the  homogeneous  (average)  velocity  gradient  of  the  mi¬ 
crostructure  and  hij(x)  is  the  velocity  gradient  fluctuation  at  point  x. 

To  model  the  deformation  of  the  microstructure,  an  irregular  material  grid 
is  expected.  However,  this  irregular  grid  in  the  real  space  results  in  difficulties 
on  conducting  Fourier  transform  in  the  next  time  step  in  a  time-dependent 
simulation,  because  that  fast  Fourier  transform  requires  a  regular  grid.  To 
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Figure  2:  (a)  The  initial  randomly  placed  Voronoi  cell  centroids  with  corresponding  circular 
grain  envelopes,  (b)  The  relaxed  centroids  which  are  used  for  Voronoi  tessellation.  The 
microstructure  is  assumed  to  be  periodic.  Numbers  in  figures  are  grain  IDs. 

resolve  this  complexity,  a  strategy  of  using  two  grids,  a  regular  computation 
grid  and  an  irregular  material  configuration  grid,  proposed  in  [43]  based  on 
the  Particle- In- Cell  (PIC)  method  [48,  49]  is  introduced.  The  computation 
grid  is  used  for  applying  the  fast  Fourier  transform  method  to  evaluate  the 
strain  related  helds.  It  is  a  regular  grid  but  not  necessary  rectangular  as 
required  by  Fourier  transform.  The  material  grid,  on  the  other  hand,  is 
attached  to  material  particles,  on  which  constitutive  relations  are  carried 
out. 

Each  grid  carries  its  own  set  of  unknowns.  Information  needs  to  be  trans¬ 
ferred  back  and  forth  between  the  two  grids  during  computation.  At  the  be¬ 
ginning  of  each  time  step,  the  initial  guess  of  local  stress  and  polarization  are 
computed  on  the  material  grid  given  the  initial  strain  rate.  The  polarization 
held  is  transferred  to  the  computation  grid,  on  which  fast  Fourier  transform 
is  performed.  The  updated  velocity  gradient  is  then  transferred  back  to  the 
material  grid  so  that  the  stress  related  helds  can  be  updated  using  the  con¬ 
stitutive  model.  At  the  end  of  the  time  step,  the  computation  and  material 
grids  are  deformed.  The  material  particles  move  according  to  local  velocity 
(Eq.  (58))  and  the  regular  computation  grid  evolves  with  the  average  velocity 
gradient: 


x(X)  =  (I  +  LAt)  X.  (59) 

A  schematic  showing  the  operation  of  the  multi-grid  strategy  is  depicted  in 

Fig.  3. 
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Figure  3:  A  schematic  description  of  the  multi-grid  CEPFFT  strategy.  The  constitutive 
model  is  applied  only  on  the  material  grid,  whereas  fast  Fourier  transform  operates  on  the 
computation  grid. 


The  information  transfer  from  regular  computation  grid  to  material  points 
is  performed  through  an  interpolation  operation.  A  family  of  shape  functions 
Af“(x),a  =  are  introduced.  Since  the  regular  computation  grid 

forms  quadrilateral  (in  2D)  or  brick  (in  3D)  elements,  a  natural  choice  of 
are  the  classical  finite  element  shape  functions.  Given  nodal  values  /c(x“) 
of  a  function  /(x)  on  the  computation  grid  (denoted  by  subscript  c),  the 
interpolated  value  at  any  material  point  Xp  in  the  microstructure  (denoted 
by  subscript  p)  is  given  as 


A 

fpi^p)  =  5^^“(xp)/c(x;?).  (60) 

a 

The  interpolation  is  local  since  hnite  element  shape  functions  are  adopted. 
Only  nodes  belonging  to  the  element  that  contains  particle  Xp  are  involved. 

The  inverse  transfer  from  the  material  grid  to  the  computation  grid  is 
performed  using  the  same  interpolation  functions  [43] .  We  assume  all  mate¬ 
rial  particles  have  the  same  mass  rup.  The  mass  of  a  computational  node  is 
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defined  as 


K 

mc(xe)  =  (61) 

k 

The  sum  is  over  all  material  particles  that  are  contained  in  the  elements 
that  share  a  common  computation  node  x^.  N°'  is  the  corresponding  shape 
function  that  connects  particle  x^  to  node  x^.  The  function  value  fc  of  /  at 
any  computation  node  x^  is  therefore: 


/c(Xc)  = 


mc(xc 


K 


;)fp 


(62) 


A  2D  illustration  of  initial  and  deformed  computation  and  material  grids 
is  shown  in  Fig.  4,  where  big  red  spots  represent  material  particles  and  small 
black  spots  are  the  nodes  of  the  computation  grid.  The  deformed  computa¬ 
tion  grid  remains  regular  so  that  fast  Fourier  transform  can  be  conducted, 
while  material  particles  move  heterogeneously.  In  most  PIC  studies,  the  ma¬ 
terial  grid  is  hner  than  the  computational  one.  In  the  current  work,  however, 
we  employ  the  same  resolution  of  the  two  grids. 


Figure  4:  (a)  Initial  computation  and  material  grids.  (2)  Deformed  computation  and 
material  grids.  The  red  dots  denote  material  particles  and  the  black  dots  denote  nodes 
of  the  computation  grid.  The  solid  black  lines  show  the  connection  between  one  material 
particle  {A)  and  computation  nodes  and  the  red  dashed  lines  are  the  connections  between 
one  computation  node  (3)  and  surrounding  material  particles. 


This  multi-grid  strategy  requires  extra  computation  time.  In  the  litera¬ 
ture,  one  regular  grid  is  used  and  the  material  grid  is  assumed  to  coincide 
with  the  computation  grid  at  all  times  [38,  40].  In  the  current  work,  we  will 
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conduct  simulations  primarily  using  this  single-grid  simplification.  Compar¬ 
ison  with  the  multi-grid  results  are  demonstrated.  It  will  be  shown  that  for 
the  examples  considered,  the  mechanical  responses  computed  from  the  multi- 
or  single-grid  approaches  do  not  vary  significantly  although  the  deformed  mi¬ 
crostructures  have  distinct  geometry. 

The  orientations  of  crystals  also  evolve  with  deformation.  Upon  conver¬ 
gence,  the  local  crystallographic  lattice  rotations  can  be  updated  by  the  spin 
rate  calculated  by 

(l»(x)  =  Q  +  (h(x)  —  (h®^*^(x),  (63) 

where  O  =  antisymiV^)  is  the  average  spin  rate  over  the  microstructure 
domain.  <-i;(x)  =  th'^(x)  -|-a;^(x)  is  the  spin  rate  fluctuation  induced  by  elastic 
and  plastic  rotation.  The  last  term  that  is  subtracted  is  the  rotation  rate  due 
to  plastic  shear  (slip)  that  does  not  distort  the  crystal  lattice.  It  is  calculated 
by 

Ns 

^*iiP(x)  =  .^U),  (64) 

Q: 

where  is  the  anti-symmetric  Schmid  tensor  (/3*'“^  =  antisym{S^°‘^)  = 
®  ®  8*^“)),  and  are  slip  direction  and  normal  to  the 

slip  plane  of  the  a-th  slip  system,  respectively). 

5.  Crystal  Plasticity  Finite  Element  Method 

We  will  compare  the  CEPFFT  results  with  conventional  full-field  crystal 
plasticity  finite  element  method  adopting  a  total  Lagrangian  scheme  [27] 
developed  in  our  previous  work  [29].  The  reference  configuration  of  each 
time  step  is  the  original  undeformed  microstructure.  The  total  deformation 
gradient  F  is  multiplicatively  decomposed  into  an  elastic  part  F®  and  a  plastic 
part  F^. 

F  =  F^F^*,  (65) 

where  plastic  deformation  is  volume  conserved:  |Fp|  =  1. 

The  velocity  gradient  of  plastic  deformation  is  determined  by  slip  rates 
on  all  slip  systems: 

LP  =  pppp-i  =  0  n(“).  (66) 
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The  plastic  strain  rate  and  spin  rate  are  calculated  by  =  symiJJ’)  and 
ijj^  =  antisym(LP) ,  respectively.  The  same  flow  rule  (e.g.  Eq.  (39))  and 
hardening  law  (e.g.  Eq.  (40))  as  in  the  CEPFFT  model  are  employed  for 
estimating  shearing  rates  7^“^ 

The  orientation  computed  at  Gauss  points  is  updated  using  the  elastic 
deformation  gradient  as: 

n(a)  _  (67) 

where  R®  is  the  rotation  part  of  F®,  and  Sq“^  and  iiq^^  are  slip  direction  and 
plane  normal  of  the  a-th  slip  system  of  the  undeformed  grain. 

6.  Numerical  examples 

In  this  section,  numerical  examples  conducted  using  the  fast  Fourier  trans¬ 
form  based  approach  are  presented.  Comparison  between  different  formula¬ 
tions  (e.g.  visco-plasticity  vs.  elasto-plasticity)  and  between  different  meth¬ 
ods  (finite  element  method  vs.  fast  Fourier  transform  based  methods)  are 
conducted  to  validate  the  current  developments.  In  addition,  the  basic  FFT 
methodology  presented  earlier  is  compared  with  the  augmented  Lagrangian 
approach  introduced  in  [40,  33].  Mechanical  response  and  fatigue  proper¬ 
ties  measured  by  strain  based  fatigue  indicator  parameters  [50]  of  INIOO 
microstructures  are  studied  using  the  novel  CEPFFT  method.  The  compu¬ 
tational  efficiency  of  CEPFFT  and  crystal  plasticity  finite  element  method 
are  presented  to  show  the  merit  of  the  CEPFFT  method.  The  use  of  the 
multi-grid  strategy  is  also  discussed. 

6.1.  Basic  formulation  versus  the  augmented  Lagrangian  formulation 

We  will  start  with  a  benchmark  plane  strain  example  of  3D  polycrys¬ 
talline  microstructure  simulated  using  the  crystal  visco-plasticity  fast  Fourier 
transform  method  (elastic  response  is  neglected).  The  crystal  visco-plastic 
constitutive  model  is  implemented  in  the  basic  framework  highlighted  ear¬ 
lier  in  this  paper  as  well  as  using  the  so  called  augmented  Lagrangian  for¬ 
mulation  [40,  33].  The  single-grid  strategy  is  adopted  so  that  deformation 
fluctuations  are  neglected. 

FCC  aluminum  is  considered.  A  cubic  polycrystalline  microstructure 
composed  of  64  grains  is  generated  using  the  Voronoi  tessellation  scheme  in 
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an  Imm^  domain.  The  microstrnctnre  is  discretized  by  16  x  16  x  16  eqnally 
spaced  voxels.  The  macroscale  velocity  gradient  is 


L  =  VV  = 


0.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

0.0 

-1.0 

X  10“^(s“^). 


(68) 


The  rate-dependent  flow  rule  (Eq.  (39))  is  used  with  70  =  ls“^  and 
m  =  0.1.  The  hardening  model  given  in  Eq.  (40)  is  adopted  to  update 
slip  resistances.  The  parameters  in  the  hardening  law  are  selected  according 
to  [40]:  kq  =  47.0MPa,  Ki  =  SQ.OMPa,  6*0  =  SSO.OMpa,  and  9i  =  IQ.OMPa. 
An  arbitrary  random  texture  is  assigned  to  the  microstrnctnre.  The  initial 
microstrnctnre  conflguration  and  pole  flgures  showing  the  random  orientation 
distribution  are  shown  in  Fig.  5. 


Figure  5:  (a)  The  image  representation  of  a  3D  polycrystalline  microstructure  containing 
64  grains,  (b)  Pole  figures  of  the  microstructure  with  randomly  assigned  orientations. 


This  example  is  used  to  validate  the  implementation  of  the  general  EFT 
framework  by  providing  a  comparison  with  the  alternative  augmented  La- 
grangian  implementation.  After  the  thickness  of  the  microstrnctnre  reduces 
by  50%,  the  deformed  microstrnctnre  and  its  stress  and  strain  fields  are  plot¬ 
ted  in  Fig.  6.  For  both  implementations,  the  stress  distribution  over  the 
microstrnctnre  is  consistent  with  the  grain  geometry.  The  local  mechani¬ 
cal  responses  of  the  two  simulations  are  very  close.  Both  intergranular  and 
intragranular  heterogeneities  of  the  stress  and  strain  rate  fields  are  captured. 

The  volume-averaged  effective  stress-strain  responses  of  the  entire  mi¬ 
crostructure  computed  by  the  two  different  algorithms  are  shown  in  Fig.  7(a). 
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(a)  (b) 


Figure  6:  Contour  plots  of  plane  strain  deformed  microstructures  evaluated  by  different 
algorithms.  The  top  layer  is  the  equivalent  (plastic)  strain  field,  and  the  bottom  layer 
is  the  equivalent  stress  field,  (a)  Crystal  visco-plasticity  fast  Fourier  transform  approach 
implemented  in  the  basic  formulation  (b)  Crystal  visco-plasticity  fast  Fourier  transform 
approach  implemented  in  the  augmented  Lagrangian  formulation. 


The  two  curves  overlap.  Pole  figures  showing  the  orientation  distribution  of 
deformed  microstructure  are  depicted  in  Fig.  7(b),  from  which  we  observe 
the  typical  plane  strain  deformation  texture  pattern  for  both  cases. 

From  the  above  comparison  of  the  local  and  effective  mechanical  re¬ 
sponses,  we  conclude  that  consistent  results  are  obtained  from  the  two  al¬ 
gorithms.  An  error  analysis  study  for  this  problem  also  has  revealed  that 
the  equilibrium  error  has  similar  convergence  for  both  approaches.  This  fact 
proves  that  there  is  no  need  for  an  augmented  Lagrangian  implementation  of 
the  FFT  algorithm  and  that  the  basic  formulation  has  sufficiently  accurate 
performance  for  solving  polycrystal  plasticity  problems  [51]. 

It  is  also  worth  mentioning  that  the  equilibrium  error  is  mostly  deter¬ 
mined  by  the  resolution  of  the  microstructure.  For  high  resolution,  the  equi¬ 
librium  condition  is  fulfilled  with  smaller  error  (see  Fig.  8).  When  the  number 
of  pixels  per  side  is  doubled,  the  equilibrium  error  is  approximately  halved. 
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Figure  7:  (a)  The  macroscopic  effective  stress-strain  responses  computed  by  the  basic  and 
augmented  Lagrangian  crystal  visco-plasticity  FFT  algorithms,  (b)  Pole  figures  of  the 
deformed  microstructure  texture. 


6.2.  Crystal  elasto -plastic  FFT  simulations  for  poly  crystalline  microstruc¬ 
tures 

We  consider  here  the  same  example  as  in  Section  6.1  nsing  the  CEPFFT 
method.  Both  the  microstrnctnre  configuration  and  material  parameters 
are  identical  to  the  previous  example.  The  elastic  constants  in  the  elasto- 
plastic  model  are  chosen  to  be  Cn  =  110  x  lO^MPa,  C12  =  59  x  lO^MPa, 
C44  =  26  X  lO^MPa.  The  CEPFFT  results  are  compared  with  the  pure 
visco-plastic  computation  as  well  as  the  results  obtained  from  the  crystal 
plasticity  finite  element  method.  Both  multi-grid  and  single-grid  strategies 
are  also  adopted  and  compared. 

In  the  finite  element  simulation,  homogeneous  boundary  condition  is  ap¬ 
plied  to  the  microstrnctnre  to  drive  its  deformation  while  the  boundary  con¬ 
ditions  of  FFT-based  simulations  are  periodic.  The  homogeneous  boundary 
condition  enforces  all  boundary  nodes  to  have  the  same  deformation/velocity 
gradient  (e.g.  Eq.  (68)  for  the  current  plane  strain  problem),  but  heteroge¬ 
neous  nodal  response  inside  the  microstrnctnre  is  allowed.  This  homogeneous 
boundary  condition  will  result  in  different  local  mechanical  responses  on  the 
microstrnctnre  from  those  obtained  using  the  periodic  boundary  condition, 
but  the  macroscopic  effective  response  should  be  comparable.  The  FEM 
simulation  is  conducted  using  our  in-house  solver  extended  based  on  [29]. 

We  first  adopt  the  main  CEPFFT  proposed  in  Section  2.4  with  single-grid 
strategy  to  perform  simulations.  The  total  strain,  plastic  strain,  and  stress 
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Figure  8:  Equilibrium  error  as  a  function  of  resolution  (number  of  pixels  per  side).  The 
equilibrium  error  is  evaluated  when  the  convergence  error  reaches  below  10“^. 

fields  of  the  deformed  microstructure  after  50%  thickness  reduction  computed 
by  different  models  are  plotted  in  Fig.  9.  It  is  seen  that  CEPFFT  gives  very 
consistent  prediction  to  both  strain  and  stress  helds  with  the  pure  visco¬ 
plastic  results  obtained  in  Section  6.1.  The  magnitude  of  the  plastic  strain 
field  prediction  by  CEPFFT  is  slightly  smaller  than  the  one  predicted  by 
visco-plastic  approach.  On  the  other  hand,  the  fields  predicted  by  the  FFT- 
based  simulations  show  differences  (especially  for  strain  fields)  from  those 
based  on  crystal  plasticity  finite  element  simulations.  The  major  causes  of  the 
differences  between  the  FFT-based  methods  and  the  FEM  approach  include: 
(1)  the  different  boundary  conditions  (periodic  by  the  FFT-based  methods 
and  homogeneous  by  FEM),  (2)  the  grain  in  FFT-based  simulation  is  assigned 
to  the  digital  model  in  the  unit  of  “point”  while  it  is  assigned  to  the  finite 
element  model  in  the  unit  of  “element”  (each  element  contains  8  points), 
(3)  the  fast  Fourier  transform  results  are  computed  directly  on  the  voxel 
points  while  the  finite  element  results  are  computed  on  integration  points 
and  extrapolated  to  nodes  using  a  least  squares  method,  and  (4)  algorithmic 
differences.  However,  we  notice  that  the  stress  fields  exhibit  similar  patterns. 
This  is  because  the  mean  strains  predicted  by  the  two  methods  are  about 
the  same  and  much  larger  in  value  than  the  strain  fluctuations.  One  thing  is 
also  worth  mentioning  is  that  the  spatial  variation  of  strain  helds  predicted 
by  crystal  plasticity  hnite  element  method  is  actually  on  the  same  level  with 
CEPFFT.  In  the  contour  plots,  the  values  shown  in  the  hnite  element  results 
is  adjusted  to  be  smaller  than  in  CEPFFT  in  order  to  show  the  stain  variation 
on  boundary  nodes,  where  homogeneous  deformation  is  applied. 
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Figure  9:  Contour  plots  of  plane  strain  deformed  microstructures  evaluated  by  different 
methods.  The  first  row  is  the  equivalent  total  strain  field,  the  second  row  is  the  equivalent 
plastic  strain  field,  and  the  bottom  row  is  the  equivalent  stress  field,  (a)  Crystal  visco¬ 
plasticity  fast  Fourier  transform  method.  The  total  strain  and  plastic  strain  are  identical 
here  since  the  elastic  response  is  ignored  in  this  model,  (b)  Crystal  elasto-plasticity  fast 
Fourier  transform  method  (CEPFFT)  (c)  Crystal  plasticity  finite  element  method. 


The  macroscopic  effective  stress-strain  responses  of  the  entire  microstruc- 
tnre  are  compared  in  Fig.  10.  The  elastic  response  is  successfnlly  captnred 
by  the  CEPFFT  model  and  is  comparable  to  the  finite  element  prediction.  It 
is  observed  that  CEPFFT  gives  close  prediction  to  the  macroscopic  response 
with  finite  element  approach  even  thongh  the  bonndary  conditions  for  the 
two  approaches  are  different. 

The  textnres  predicted  by  the  three  models  are  plotted  using  pole  figures 
in  Fig.  11.  It  is  observed  that  all  the  three  sets  of  pole  figures  predicted  are 
very  similar.  The  typical  plane  strain  pattern  is  obtained. 

To  check  the  convergence  of  the  results  with  respect  to  increasing  reso¬ 
lution,  we  performed  the  same  simulation  on  the  same  microstructure  but 
discretized  by  32  x  32  x  32  voxels  using  CEPFFT.  The  local  and  effective 
mechanical  responses  are  compared  in  Figs.  12  and  13,  respectively.  It  is  ob¬ 
served  that  the  two  sets  of  results  are  consistent.  Plotted  in  the  two  hgures 
are  also  the  predictions  using  the  modihed  CEPFFT  approach  introduced  in 
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Figure  10:  The  macroscopic  effective  stress-strain  responses  of  plane  strain  deformed  mi¬ 
crostructures  predicted  by  different  models,  (a)  Effective  stress-total  strain  responses  by 
CEPFFT  and  crystal  plasticity  FEM;  (b)  The  effective  stress-plastic  strain  responses  by 
the  three  methods.  Note  that  here  CVPEET  denotes  crystal  visco-plasticity  fast  Fourier 
method,  and  CPFEM  refers  to  crystal  plasticity  finite  element  method. 


Remark  2,  where  stress  is  linked  to  total  strain  rate  through  an  elasto-plastic 
modulus.  We  observe  that  the  results  from  this  modihed  implementation  are 
close  to  the  main  formulation  of  this  paper  (that  computes  elastic  and  plastic 
fluctuations  separately) ,  although  with  milder  spatial  variation,  implying  its 
good  performance  for  this  problem. 

We  also  studied  the  equilibrium  error  of  the  main  CEPFFT  approach 
for  microstructures  with  different  resolution.  The  improvement  of  error  with 
rehning  the  image  is  seen  from  Fig.  14(a),  which  is  similar  to  the  visco¬ 
plasticity  case.  The  convergence  of  the  CEPFFT  model  as  a  function  of 
iteration  number  is  shown  in  Fig.  14(b).  We  observe  a  fast  convergence  rate 
for  all  tests. 

In  order  to  check  the  effect  of  heterogeneous  deformation  on  mechanical 
responses,  we  next  repeat  the  above  simulations  with  the  multi-grid  strat¬ 
egy.  implemented  with  the  CEPFFT.  The  deformed  material  grid  obtained 
by  the  multi-grid  strategy  and  the  homogeneous  approximation  are  plotted 
in  Fig.  15.  Contour  plots  of  local  mechanical  responses  are  shown  in  Fig.  16. 
Comparing  with  the  single-grid  results,  we  hnd  that  the  predicted  mechanical 
responses  do  not  vary  much,  although  the  deformed  microstructures  becomes 
irregular.  The  effective  stress-strain  curve  and  grain  orientation  distribution 
predicted  by  the  multi-grid  strategy  are  almost  identical  with  the  correspond- 
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Figure  11:  Crystallographic  textures,  represented  in  pole  figures,  of  plane  strain  de¬ 
formed  microstructures  predicted  by  three  models,  (a)  Crystal  visco-plasticity  fast  Fourier 
transform  method  (CVPFFT)  (b)  Crystal  elasto-plasticity  fast  Fourier  transform  method 
(CEPFFT)  (c)  Crystal  plasticity  finite  element  method  (CPFEM). 


ing  results  obtained  with  the  single-grid  prediction  and  are  not  repeated  here. 

We  next  examine  an  example  of  volume-conserved  compression  using 
three  different  models.  The  same  polycrystalline  microstructure  contain¬ 
ing  64  grains  and  material  parameters  are  used,  while  the  imposed  velocity 
gradient  becomes 


L  =  VV 


0.5  0.0  0.0 
0.0  0.5  0.0 
0.0  0.0  -1.0 


X  10“®(s“‘). 


(69) 


The  single-grid  strategy  along  with  the  main  CEPFFT  formulation  is 
firstly  adopted.  After  the  compression  in  z  direction  reaches  50%,  the  total 
strain,  plastic  strain,  and  stress  helds  are  plotted  in  Fig.  17.  The  local  me¬ 
chanical  helds  obtained  from  visco-plasticity  and  elasto-plasticity  approaches 
agree  very  well.  We  also  observe  that  the  hnite  element  stress  and  strain  held 
results  show  very  similar  pattern  with  the  FFT-based  results  despite  the  dif¬ 
ferent  kinematic  conditions  imposed  on  the  FEM  and  FFT  methodologies. 


The  macroscopic  effective  stress-strain  responses  of  the  entire  microstruc¬ 
ture  are  compared  in  Fig.  18.  The  elastic  response  is  also  successfully  cap¬ 
tured  by  the  CEPFFT  model  and  is  comparable  with  the  hnite  element 
prediction. 

The  crystallographic  textures  described  by  pole  hgures  are  also  captured 
(Fig.  19).  We  can  see  once  more  that  the  diherent  models  give  very  close 
prediction  on  the  texture  evolution  of  the  volume-conserved  compression. 

To  check  the  convergence  of  the  results  with  respect  to  increasing  resolu¬ 
tion  for  volume-conserved  compression  deformation,  we  also  performed  the 
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Figure  12:  Contour  plots  of  plane  strain  deformed  microstructures  with  different  resolution 
and  methods.  The  first  row  is  the  equivalent  total  strain  field,  the  second  row  is  the 
equivalent  plastic  strain  field,  and  the  bottom  row  is  the  equivalent  stress  field,  (a)  The 
main  CEPFFT  method  using  16  x  16  x  16-voxel  microstructure,  (b)  The  main  CEPFFT 
method  using  32  x  32  x  32-voxel  microstructure,  (c)  The  modified  CEPEET  (Remark  2) 
using  16  x  16  X  16-voxel  microstructure. 


same  simulation  on  the  same  microstructure  but  discretized  by  32  x  32  x  32 
voxels.  The  local  and  effective  mechanical  responses  are  compared  in  Figs.  20 
and  21,  respectively.  The  results  computed  using  the  modified  CEPFFT  for¬ 
mulation  discussed  in  Remark  2  are  also  demonstrated.  Both  the  local  and 
effective  responses  are  close  for  the  two  implementations  of  the  CEPFFT 
method. 

Comparing  the  multi-grid  with  the  single-grid  results,  we  note  that  the 
irregular  deformation  is  captured  (Fig.  22),  while  the  mechanical  responses  do 
not  exhibit  obvious  change  (Fig.  23).  Almost  identical  effective  stress-strain 
response  and  texture  evolution  with  the  single-grid  prediction  are  observed 
and  as  such  these  results  are  not  shown  here. 

6.3.  Investigation  of  fatigue  indicator  parameters  of  IN  100 

In  this  subsection,  we  adapt  the  CEPFFT  method  to  study  fatigue  prop¬ 
erties  of  Ni-based  superalloys  using  the  flow  rule  and  hardening  law  for  the 
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Figure  13:  (a)  The  macroscopic  effective  stress-total  strain  responses  by  16  x  16  x  16 
microstructure  and  32  x  32  x  32  microstructure  obtained  using  different  formulations, 
(b)  Crystallographic  textures  represented  in  pole  figures.  Main  CEPFFT  refers  to  the 
main  crystal  elasto-plasticity  FFT  method  and  Modified  CEPFFT  refers  to  the  crys¬ 
tal  elasto-plasticity  implementation  using  a  homogeneous  elasto-plastic  medium  approach 
(Remark  2). 


INIOO  model  introduced  in  Section  3.1.2.  The  same  microstructure  config¬ 
uration  as  in  the  previous  examples  is  used  and  the  initial  texture  is  ran¬ 
dom.  The  same  parameters  of  the  constitutive  equations  as  listed  in  [30] 
are  used.  The  volume  fractions  and  sizes  of  7'  precipitates  are  given  by 
fpi  =  0,  fp2  =  0.42,  d2  =  lOSnm,  fp^  =  0.11,  =  7nm.  The  microstructure 

is  subjected  to  a  3-loop  cyclic  loading  (tension  and  compression  along  the 
^-direction) .  The  stress-strain  response  in  the  ^-direction  during  the  3  load¬ 
ing  loops  is  plotted  in  Fig.  24  with  comparison  to  an  FEM  simulation  that 
adopts  the  same  constitutive  model.  The  CEPFFT  models  implemented  in 
both  the  main  formulation  and  the  modified  formulation  (in  Remark  2)  are 
tested.  The  single-grid  update  strategy  is  adopted.  The  microstructure  in¬ 
put  to  the  CEPFFT  simulation  is  discretized  by  16  x  16  x  16  voxels  to  be 
consistent  with  the  finite  element  input  (16  x  16  x  16  cubic  elements).  Homo¬ 
geneous  boundary  condition  is  applied  to  the  finite  element  microstructure. 
The  CEPFFT  simulations  give  similar  prediction  to  the  stress-strain  “loop” 
with  the  finite  element  model.  Note  that  the  loop  predicted  by  the  main 
CEPFFT  algorithm  is  wider  than  the  loops  obtained  from  the  FEM  or  the 
modified  CEPFFT  algorithm. 
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Figure  14:  (a)  Equilibrium  error  as  a  function  of  resolution  (number  of  pixels  per  side). 
The  equilibrium  error  is  evaluated  using  the  basic  formulation  when  the  convergence  error 
reaches  below  10“^.  (b)  Convergence  error  versus  number  of  iterations. 


(a)  (b) 


Figure  15:  (a)  Deformed  microstructure  predicted  by  multi-grid  CEPFFT.  (c)  Deformed 
microstructure  predicted  by  single-grid  CEPFFT. 


We  next  compute  strain  based  fatigue  indicator  parameters  (FIPs)  related 
to  small  crack  formation  and  early  growth  [50].  The  three  fatigue  indica¬ 
tor  parameters  of  interest  are  the  cumulative  plastic  strain  per  cycle  {Pcyc), 
which  correlates  to  the  crack  incubation  life;  the  Fatemi-Socie  parameter 
(Pps),  which  relates  to  the  small  crack  growth;  and  the  maximum  range  of 
cyclic  plastic  shear  strain  parameter  {Pmps)-  The  dehnitions  of  these  fatigue 
indicator  parameters  can  be  found  in  [44].  Contour  plots  of  fatigue  indicator 
parameters  over  the  microstructure  are  plotted  in  Fig.  25.  All  fatigue  indica¬ 
tor  parameters  are  computed  for  the  last  (3rd)  loading  loop.  We  observe  that 
the  fatigue  indicator  parameters  contour  plots  demonstrate  similar  patterns 
in  CEPFFT  and  crystal  plasticity  hnite  element  simulations.  The  results 
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(a)  (b)  (c) 


Figure  16:  Contour  plots  of  plane  strain  deformed  microstructures  computed  using  dif¬ 
ferent  strategies.  The  first  row  shows  results  obtained  from  the  multi-grid  CEPFFT,  and 
the  bottom  row  shows  results  obtained  from  the  single-grid  CEPFFT.  (a)  Equivalent  total 
strain  (b)  Equivalent  plastic  strain  (c)  Equivalent  stress. 


predicted  by  the  main  CEPFFT  approach  show  more  heterogeneity  and  are 
closer  to  hnite  element  results  than  the  modihed  FFT  model  (Remark  2). 

The  convergence  of  the  CEPFFT  simulation  with  respect  to  resolution 
and  the  multi-grid  effect  are  also  tested  for  this  fatigue  problem.  The  fa¬ 
tigue  indicator  parameters  are  plotted  in  Fig.  26.  It  is  observed  that  the 
deformation  heterogeneity  does  not  affect  significantly  the  fatigue  indicator 
parameters  helds  and  the  distortion  of  the  microstructure  is  small  since  the 
total  strain  is  small  for  this  problem.  The  hner  microstructure  (32  x  32  x  32 
voxels)  gives  consistent  prediction  to  local  fatigue  indicator  parameters  with 
the  coarse  grid  prediction,  while  larger  heterogeneity  is  captured. 

Grain  level  fatigue  indicator  parameters,  namely  the  maximum  and  av¬ 
erage  fatigue  indicator  parameters  of  voxels  within  individual  grains  [44]  are 
also  extracted.  Their  distributions  are  shown  in  Fig.  27.  The  distributions 
computed  using  the  main  CEPFFT  based  on  the  16  x  16  x  16-voxel  mi¬ 
crostructure  are  shown  in  (a).  The  ones  computed  using  the  main  CEPFFT 
on  the  32  X  32  X  32- voxel  microstructure  are  plotted  in  (b).  The  results  of  the 
modihed  implementation  of  CEPFFT  (Remark  2)  are  depicted  in  (c),  and 
the  crystal  plasticity  hnite  element  results  are  given  in  (d).  These  results  are 
obtained  using  the  single-grid  strategy.  In  the  hgure,  the  FIP  with  “max” 
prehx  means  the  maximum  value  of  the  FIP  of  voxels  within  an  individual 
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Figure  17:  Contour  plots  of  compressed  microstructures  evaluated  by  different  methods. 
The  first  row  is  equivalent  total  strain  field,  the  second  row  is  equivalent  plastic  strain 
field,  and  the  bottom  row  is  equivalent  stress  field,  (a)  Crystal  visco-plasticity  fast  Fourier 
transform  method  (CVPFFT).  The  total  strain  and  plastic  strain  are  identical  here  since 
the  elastic  response  is  ignored  in  this  model,  (b)  Crystal  elasto-plasticity  fast  Fourier 
transform  method  (CEPFFT).  (c)  Crystal  plasticity  finite  element  method  (CPFEM). 


grain.  The  one  with  “ave”  prehx  indicates  that  it  is  the  average  FIP  over 
voxels  within  an  individual  grain.  The  distributions  are  normalized  by  the 
maximum  value  of  each  FIP  over  all  grains  of  the  microstructure.  It  is  seen 
that  the  distributions  of  grain  level  fatigue  indicator  parameters  predicted  by 
the  main  CEPFFT  on  the  hner  microstructure  are  closer  to  the  hnite  element 
results.  The  ones  predicted  by  the  main  CEPFFT  on  coarse  microstructure 
are  also  close  to  the  hne  microstructure  estimations.  On  the  other  hand,  the 
maximum  and  average  fatigue  indicator  parameters  distributions  predicted 
by  the  modihed  CEPFFT  are  not  very  distinguishable,  implying  that  the 
intragranular  heterogeneity  of  fatigue  indicator  parameters  predicted  by  the 
modihed  CEPFFT  formulation  is  weak.  This  is  consistent  with  the  contour 
plots  of  the  FIP  helds.  A  more  sophisticated  elasto-plastic  modulus  may 
resolve  this  problem. 
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Figure  18:  The  macroscopic  effective  stress-strain  responses  of  compressed  microstructures 
predicted  by  different  models,  (a)  Effective  stress-total  strain  responses  by  CEPFFT 
and  crystal  plasticity  finite  element  method  (CPFEM);  (b)  Effective  stress-plastic  strain 
responses  by  all  three  methods.  CVPFFT  refers  to  crystal  visco-plasticity  fast  Fourier 
method. 


6.4-  Computation  Efficiency 

The  computational  efficiency  of  CEPFFT  method  is  important  since  our 
goal  is  to  develop  a  high  efficient  full-held  simulator  that  can  be  adopted  in 
stochastic  simulations  as  the  deterministic  solver  to  further  study  the  proba¬ 
bilistic  nature  of  material  properties.  As  mentioned  before,  with  the  employ¬ 
ment  of  fast  Fourier  transform  to  solve  for  local  mechanical  responses,  the 
cumbersome  matrix  inversion  in  hnite  element  simulation  is  circumvented, 
which  leads  to  great  improvement  of  the  computation  speed.  The  computa¬ 
tion  times  of  microstructures  subjected  to  plane  strain  deformation  to  0.02 
strain  are  shown  in  Fig.  28(a)  and  tabulated  in  Table  1.  All  tests  were  run  on 
the  Teragrid  TACC  Lonestar  Linux  Cluster.  Each  computation  node  contains 
two  Xeon  Intel  Hexa-Core  64-bit  Westmere  processors,  the  core  frequency  of 
which  is  3.33GHz.  Single-grid  morphology  evolution  strategy  and  main  for¬ 
mulation  are  adopted  by  all  CEPFFT  tests.  We  hrst  ran  examples  on  one 
single  processor,  respectively.  The  CEPFFT  simulation  of  a  microstructure 
discretized  by  16^  voxels  only  takes  less  than  3  minutes,  while  the  hnite  ele¬ 
ment  simulation  of  a  microstructure  containing  16^  elements  takes  13  hours. 
With  an  increased  resolution  of  32^  voxels,  the  computation  time  of  CEPFFT 
increases  to  about  22  minutes,  which  is  still  signihcantly  smaller  than  the  h- 
nite  element  simulation.  In  order  to  further  reduce  the  computation  time. 
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Figure  19:  Crystallographic  textures,  represented  in  pole  figures,  of  compressed  mi¬ 
crostructures  predicted  by  three  models,  (a)  Crystal  visco-plasticity  fast  Fourier  trans¬ 
form  method  (CVPFFT).  (b)  Crystal  elasto-plasticity  fast  Fourier  transform  method 
(CEPFFT).  (c)  Crystal  plasticity  finite  element  method  (CPFEM). 


Table  1:  Computation  times  for  microstructures  under  plane  strain  and  cyclic  deforma¬ 
tions  simulated  using  different  methods.  CPFEM  refers  to  crystal  plasticity  finite  element 
method. _ 


1-processor, 
CEPFFT, 
16^- voxel 

60-processor, 
CEPFFT, 
16^- voxel 

1-processor, 

CEPFFT, 

32^-voxel 

240-processor, 

CEPFFT, 

32^-voxel 

1-processor, 

CPFEM, 

16^-element 

240-processor, 

CPFEM, 

16^-elenient 

Plane  strain 

164s 

5s 

1309s 

17s 

46389s 

560s 

Cyclic  loading 

554s 

16s 

4385s 

57s 

46873s 

522s 

we  also  parallelized  our  code.  The  16^-voxel  CEPFFT  simulation  only  takes 
5  seconds  to  finish  if  run  on  60  processors.  The  32^-voxel  simulation  running 
on  240  processors  takes  17  seconds  to  finish.  On  the  other  hand,  the  finite 
element  simulation  of  16^-element  microstructure  takes  about  10  minutes  to 
complete  when  240  processors  are  used. 

Besides  plane  strain  deformation,  we  also  tested  the  computational  ef¬ 
ficiency  of  the  microstructure  fatigue  problem.  The  computation  times  of 
one  complete  loading  loop  are  captured  and  plotted  in  Fig.  28(b).  We  ob¬ 
serve,  again,  that  the  CEPFFT  simulation  is  much  faster  than  the  crystal 
plasticity  finite  element  method.  Note  that  the  time  increment  required  by 
the  finite  element  simulation  is  smaller  than  that  required  by  the  fast  Fourier 
transform  simulation  in  order  to  reach  convergence.  However,  even  if  the  two 
approaches  use  the  same  time  increment,  the  FFT  computation  is  still  much 
faster  than  the  FEM  simulation  for  the  microstructure  resolution. 


With  the  adoption  of  fast  Fourier  transform  for  solving  the  polycrystalline 
microstructure  boundary  value  problem,  the  most  time  consuming  operation 


41 


Figure  20:  Contour  plots  of  compressed  microstructures  with  different  resolution  and 
formulation.  The  first  row  is  the  equivalent  total  strain  field,  the  second  row  is  the  equiv¬ 
alent  plastic  strain  field,  and  the  bottom  row  is  the  equivalent  stress  field,  (a)  The  main 
CEPFFT  using  16  x  16  x  16-voxel  microstructure,  (b)  The  main  CEPFFT  method  using 
32  X  32  X  32-voxel  microstructure  (c)  The  CEPFFT  method  implemented  using  a  homoge¬ 
neous  elasto-plastic  medium  approach  (Remark  2)  using  16  x  16  x  16- voxel  microstructure. 


in  the  FFT-based  simulation  is  then  the  nonlinear  constitutive  calculations 
dehned  by  the  crystal  plasticity  theory.  Another  branch  of  research  using 
Fourier  methods  in  computational  crystal  plasticity  has  been  conducted  by 
Kalidindi  and  coworkers,  who  have  introduced  Fourier  expansion  to  establish 
hypersurfaces  of  constitutive  equations  so  that  the  iterative  computation  in 
updating  nonlinear  local  properties  can  be  approximated  using  spectral  in¬ 
terpolation  [52,  53,  54,  55].  Different  from  our  strategy  of  solving  the  over¬ 
all  microstructure  governing  equations  using  Fourier  transform,  they  replace 
nonlinear  local  constitutive  equations  by  Fourier  expansion  after  calibrating 
coefficients  of  the  Fourier  series  using  a  given  database.  The  computation 
efficiency  of  a  crystal  plasticity  problem  may  be  further  improved  by  cou¬ 
pling  their  spectral  method  with  the  present  crystal  plasticity  fast  Fourier 
transform  framework. 
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Figure  21:  (a)  The  macroscopic  effective  stress-total  strain  responses  by  16  x  16  x  16 
microstructure  and  32  x  32  x  32  microstructure  obtained  using  different  formulations, 
(b)  Crystallographic  textures  represented  in  pole  figures.  Main  CEPFFT  refers  to  the 
main  crystal  elasto-plasticity  FFT  method  and  Modified  CEPFFT  refers  to  the  crys¬ 
tal  elasto-plasticity  implementation  using  a  homogeneous  elasto-plastic  medium  approach 
(Remark  2). 


7.  Conclusions 

In  this  work,  we  developed  an  efficient  FFT  fnll  field  model  to  inves¬ 
tigate  elasto-plastic  properties  of  polycrystalline  materials  by  interrogating 
image-based  realistic  microstructnres.  The  elastic  and  plastic  responses  were 
compnted  separately.  An  integrated  formulation  was  also  proposed  using  a 
particular  choice  for  the  elasto-plastic  moduli.  Predictive  capability  and  com¬ 
putational  efficiency  of  the  newly  developed  CEPFFT  method  were  presented 
using  numerical  examples  in  comparison  with  visco-plasticity  approach  and 
crystal  plasticity  finite  element  simulations.  Error  tests  were  conducted  to 
show  the  comparable  performance  of  the  basic  and  augmented  Lagrangian 
algorithms.  The  multi-grid  strategy  was  implemented  to  allow  comparison 
with  the  single-grid  simplification.  The  self-consistence  (convergence  with 
refining  discretization)  of  the  CEPFFT  method  was  shown  through  simula¬ 
tions  with  increasing  resolution  of  the  microstructure.  Fatigue  properties  of 
Ni-based  superalloy  microstructnres  described  by  fatigue  indicator  parame¬ 
ters  were  studied  using  the  CEPFFT  method.  The  main  discoveries  of  this 
work  are  as  follows: 
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(a)  (b) 


Figure  22:  (a)  Deformed  microstructure  predicted  by  multi-grid  CEPFFT.  (c)  Deformed 
microstructure  predicted  by  single-grid  CEPFFT. 

1.  The  elastic  response  of  the  microstructure  is  successfully  captured  by 
CEPFFT.  Fatigue  properties  of  nickel-based  superalloy  microstructures 
can  be  efficiently  investigated  using  CEPFFT  method.  The  computed 
results  are  comparable  with  the  those  based  on  crystal  plasticity  hnite 
element  simulation. 

2.  The  equilibrium  error  resulted  from  the  augmented  Lagrangian  algo¬ 
rithm  is  comparable  with  the  one  resulting  from  the  basic  formulation 
and  depends  on  the  resolution  of  the  input  image. 

3.  The  CEPFFT  model  provides  consistent  prediction  of  the  macroscopic 
effective  mechanical  responses  and  local  stress  held  with  the  hnite  el¬ 
ement  and  crystal  visco-plasticity  fast  Fourier  transform  simulations. 
The  crystallographic  texture  patterns  of  microstructures  under  diherent 
deformation  modes  estimated  by  the  three  methods  are  almost  identi¬ 
cal. 

4.  The  local  strain  helds  obtained  by  the  CEPFFT  implementation  agree 
very  well  with  visco-plastic  results  and  have  similar  patterns  with  the 
crystal  plasticity  hnite  element  results. 

5.  The  multi-grid  strategy  that  uses  separate  computation  and  mate¬ 
rial  grids  was  shown  to  successfully  capture  irregular  conhguration  of 
deformed  microstructures.  However,  it  does  not  result  in  signihcant 
changes  in  the  mechanical  response  for  the  tested  examples  in  compar¬ 
ison  to  the  single-grid  method. 

6.  The  required  computation  time  of  CEPFFT  is  signihcantly  less  than 
that  of  crystal  plasticity  hnite  element  method  since  the  methodology 
does  not  require  the  solution  of  any  matrix  equations  as  in  the  FEM. 
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(a)  (b)  (c) 


Figure  23:  Contour  plots  of  compressed  microstructures  computed  using  the  multi-grid 
strategy.  The  first  row  shows  results  obtained  from  the  multi-grid  CEPFFT,  and  the 
bottom  row  shows  results  obtained  from  the  single-grid  CEPFFT.  (a)  Equivalent  total 
strain,  (b)  Equivalent  plastic  strain,  (c)  Equivalent  stress. 


This  computational  efficiency  provides  significant  potential  in  integrat¬ 
ing  the  FFT  approach  with  stochastic  multiscale  materials  simulations. 

7.  Various  CEPFFT  can  be  introduced.  The  approach  highlighted  in  this 
paper  computes  separately  the  elastic  and  plastic  responses.  Care¬ 
ful  design  of  the  elasto-plastic  modulus  for  the  homogeneous  reference 
medium  can  improve  the  computational  efficiency. 
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Figure  24:  Stress-strain  responses  of  three  loops  of  cyclic  loading  computed  by  CEPFFT 
and  CPFEM.  Main  CEPFFT  refers  to  the  main  crystal  elasto-plasticity  FFT  method, 
modified  CEPFFT  refers  to  the  crystal  elasto-plasticity  implementation  using  a  homoge¬ 
neous  elasto-plastic  medium  approach  (Remark  2),  and  CPFEM  is  the  crystal  plasticity 
finite  element  method. 


Figure  25:  Contour  plots  of  fatigue  indicator  parameters  fields.  Main  CEPFFT  results 
are  placed  in  the  top  row,  results  from  the  modified  CEPFFT  formulation  (Remark  2)  are 
placed  in  the  middle  row,  and  crystal  plasticity  finite  element  results  are  located  in  the 
bottom  row.  (a)  Pcyc,  (b)  Pps,  (c)  Pmps- 
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Figure  26:  Contour  plots  of  fatigue  indicator  parameters  fields  evaluated  by  the  multi-grid 
CEPFFT  on  coarse  microstructure  (top  row),  single-grid  CEPFFT  on  coarse  microstruc¬ 
ture  (middle  row),  and  single-grid  CEPFFT  on  fine  microstructure  (bottom  row),  (a) 
PcyC7  (b)  PfS,  (c)  Pmps- 
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Figure  27:  Distribution  of  the  fatigue  indicator  parameters  among  grains  computed  by  (a) 
Main  CEPFFT  with  coarse  microstructure,  (b)  Main  CEPFFT  with  fine  microstructure, 
(c)  Modified  CEPFFT  (Remark  2)  with  coarse  microstructure  and  (d)  crystal  plasticity 
finite  element  method.  These  distributions  are  normalized  by  the  maximum  values  of  each 
FIP  over  all  grains. 


Figure  28:  Computation  times  of  simulations  using  different  methods  on  microstructures 
with  different  resolutions,  (a)  Plane  strain  deformation  to  0.02  strain,  (b)  One  complete 
loop  of  cyclic  loading  on  the  INIOO  microstructure.  The  computation  times  are  shown  in 
logarithmic  scale.  In  the  figure,  IProcCEPFFT-I6P  means  the  CEPFFT  simulation  of  a 
microstructure  discretized  by  16  x  16  x  16  voxels  using  1  processor  and  lProcCPFEM-16E 
means  the  crystal  plasticity  finite  element  simulation  of  a  microstructure  discretized  by 
16x16x16  elements  using  1  processor. 
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